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Recent microfluidic experiments revealed that large particles advected in a fluidic loop display 
long-range hydro dynamic interactions. However, the consequences of such couplings on the traffic 
dynamics in more complex networks remain poorly understood. In this letter, we focus on the 
transport of a finite number of particles in one-dimensional loop networks. By combining numerical, 
theoretical, and experimental efforts, we evidence that this collective process offers a unique example 
of Hamiltonian dynamics for hydrodynamically interacting particles. In addition, we show that the 
asymptotic trajectories are necessarily reciprocal despite the microscopic traffic rules explicitly break 
the time reversal symmetry. We exploit these two remarkable properties to account for the salient 
features of the effective three-particle interaction induced by the exploration of fluidic loops. 
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The long-range nature of the hydrodynamic interac- 
tions is responsible for fascinating collective phenomena 
in non-equilibrium suspensions, such as the velocity fluc- 
tuations of sedimenting particles [1 , and the emergence 
of coherent structures in isotropic suspensions of active 
particles However, in confined geometries, the walls 
screen exponentially the correlations of the particle veloc- 
ity fj|. Hence, no collective traffic phenomenon can oc- 
cur when dilute suspensions flow in ducts having a width 
comparable to the particle size. Nonetheless, recent mi- 
crofluidic experiments in channels including a loop, re- 
vealed a rich variety of collective dynamics, such as mul- 
tiperiodic and multistable traffic patterns fiUIO]. These 
experimental observations have been rationalized on the 
basis of two empirical rules [6' : as a particle enters a loop, 
it takes the branch in which the flow rate is maximal, 
and (ii) the particles partly obstruct the branch in which 
they flow. Consequently, the particle velocity at a node 
is a function of the particle positions in the whole loop, 
thereby inducing localized but long-range hydrodynamic 
interactions. So far, most of the research on microfluidic 
traffic flows have been dedicated to the transport through 
a single fluidic loop fed at a constant rate by a continuous 
droplet /bubble stream. 

In this letter, we investigate the dynamics of a flnite 
number of particles cruising in an extended loop-network, 
see Fig. [l] We henceforth focus on the three-body prob- 
lem. This setup is the basic building block to model the 
traffic dynamics of dilute suspensions (the case of two 
particles being trivial). Our primary idea is to consider 
the traffic through a single loop as a scattering process, 
which maps the distances \{n) — [Ai(n), A2(n)] between 
the three particles entering the loop n into a new set of 
distances A(n -I- 1) = 5[A(n)], where S is the scattering 
map. The transport through the entire network is then 
conceived as a discrete dynamical system, for which the 
loop index n stands for the time variable. From this per- 
spective, we first evidence that, remarkably, the asymp- 
totic traffic dynamics is Hamiltonian. To the best of our 




FIG. 1: Sketch of our theoretical and experimental set-up: 
three droplets are advected in a one dimensional microfluidic 
loop-network. A2 (resp Ai) is the distance between the two 
righmost (resp. leftmost) droplets. 

knowledge, this is the only system of hydrodynamically 
interacting particles, for which an Hamiltonian descrip- 
tion exists. Moreover, we show that the dynamics is 
asymptotically invariant upon time reversal symmetry 
despite the microscopic traffic rules are explicitly non- 
reciprocal. We exploit these two features to account for 
the geometrical and the dynamical properties of the scat- 
tering map S. We close this paper, by comparing our 
theoretical predictions to microfluidic experiments. A 
quantitative agreement is found without any free fltting 
parameter. 

We use a well established framework to model the traf- 
flc dynamics in a fluidic network made of a chain of N 
identical loops |7]. Precisely, it consists of four rules, 
which have proven to yield excellent agreement with the 
experiments |SVnU| : (i) The flow state of the fluid in the 
network is given by the analogous of the Kirchhoff laws, 
(ii) The particles are supposed to have a constant mo- 
bility coefficient. Therefore, we identify the fluid and 
the particle velocities, (iii) When it reaches a vertex, a 
particle takes the branch where the fluid velocity is the 
higher. Note that this empirical rule, observed on de- 
formable particles, is explicitly non-reciprocal, (iv) The 
particles partly obstruct the channels in which they jour- 
ney. Precisely, the hydrodynamic resistances, expressed 
in unit-length, are given by: L;,^(n;,^) = L;,^ -I- n^.t^p, 
where the n^^^ are the numbers of particles advected in 
the upper and in the lower branches respectively. The 
^i.t represent the branches' length, and Lp is the con- 
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FIG. 2: Numerical results obtained for: a = 10/9. A-Basins 
of attraction of the 5-map. The grey dots correspond to ini- 
tial distances yielding stationary asymptotic dynamics (zone 
C). The red dots converge to closed periodic orbits, c = 10/9. 
B- Dots: Superimposed asymptotic trajectories for c = 20/9. 
The 8 polygons correspond to the 8 trafficking scenarios. C- 
Dots: Superimposed asymptotic trajectories for c = 10/9 
(same parameters as in A). D- Close-up of the edge of one 
island, same parameters as in C. 

stant additional resistance induced by a single droplet. 
It follows that the particle velocity in the upper branch 
is: 



where v is the fluid velocity outside the loops. A sym- 
metric formula holds for the lower branch. We apply the 
above rules numerically using the event-driven algorithm 
introduced in [111 for a 1-loop network and iterate it N 
times. We record the distances Ai(n) and A2(n), irre- 
spective of the particle ordering, between the particles 
entering the loop n, see Fig. [1] The traflic dyn amies IS 
parametrized by two dimensionless numbers: the loop- 
aspect ratio a = L^/L^ and the "clogging-parameter" 
c = Lp/L^, which quantifies how much a particle hin- 
ders the flow in a given branch. In all that follows, we 
restrain ourselves to weakly asymmetric loops for which 
1 < a < 2 and a < 1 + c. Within this approximation, a 
particle entering a loop journeys through the less occu- 
pied branch, or through the lower branch if the loop is 
empty. 

The gross features of the traflic dynamics do not de- 
pend on the aspect ratio a. The phase plane (Ai,A2) is 
divided into two basins of attraction. Fig. [2]A. Starting 
from the rightmost basin, the system is quickly absorbed 
into the upper-right part of the plane, where the interpar- 
ticle distances remain constant, Figs.|2|3 and 2C. Starting 



from the leftmost basin, the system reaches a compact re- 
gion flUed with a continuous ensemble of periodic orbits 
centered on a unique fixed point. The distances oscillate 
around a constant value, and the orbits are either one 
dimensional curves, or zero dimensional (viz. returning 
repeatedly to a finite number of points). Two typical ex- 
amples are given in Figs. [2^ and|2p. These two regions 
are reached after ~ 10 loops. 

The absorbing region is defined by the inter-particle 
distance A^ax above which a particle enters the loop af- 
ter the previous one has left it. This distance is easily 
deduced from Eq. [l] A^ax = (l + ^). Above Amax, 
the particles do not interact, and the scattering-map is 
trivial whatever a and c: 5 = 1. Conversely, the geome- 
try and the topology of the closed orbits strongly depend 
on the clogging parameter. First, we shall distinguish 
two regimes from the periods of the orbits. In Fig. [3] we 
show the variations of the typical oscillation period, r, 
as a function of the clogging parameter c. There exist 
two limit values, c~ and c^, below and above which all 
the orbits share the same period. Moreover, for c < c^, 
and c > , the ID-orbits are self-similar ellipses cen- 
tered on a marginally stable fixed point, which lies on 
the line Ai = A2, see Fig. [2|3. Therefore, the iS-map is 
necessarily afhne for this range of parameters. We em- 
phasize that Ai = A2 is a symmetry axis of the ellipses 
independently of a and c. In contrast, for intermediate 
clogging parameters (c^ < c < c+), r strongly depends 
on the initial conditions, and more than one period is 
detected, see Fig. [3] In addition, at least one of the in- 
variant curves has a non-elliptical shape. However, the 
global symmetry of the phase portrait with respect to 
Ai = A2 is preserved, see Fig. [2p. Several closed orbits 
are destabilized into separatrix and island chains centered 
on stable p-periodic points. Trajectories with p = 15 are 
clearly seen in Fig. |2p. We systematically observed a hi- 
erarchy of island chains, as exemplifled in the close-up 
shown in Fig. [2p. The inner part of the largest islands 
clearly include island chains as well. They are separated 
by large chaotic regions, which also exist at the largest 
scale of the phase portrait, though they are much less 
extended. 

We close this numerical section with the flrst main re- 
sult of this letter. Remarkably, all the features of the 
phase portrait are the hallmarks of Hamiltonian map- 
pings, despite the traffic dynamics is a driven dissipative 
process. We shall note that fluid mechanics offers other 
examples of Hamiltonian descriptions for advected par- 
ticles. However, these models have so far been restricted 
to non-interacting passive tracers in bidimensional and 
incompressible fluids, for which the stream function read- 
ily provides an effective Hamiltonian [T^]. The system, 
we consider here, does not belong to this class. Both 
the loop geometry and the effective hydrodynamic cou- 
pling between the particles, make impossible the use of a 
stream function as an effective Hamiltonian. 
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To elucidate the Hamiltonian nature of the traffick- 
ing dynamics, we construct explicitly the scattering 
map S. To do so, we first note that there exist 8 
different traffic scenarios, labeled hy X — A, . . . ,H. 
These scenarios are defined by the time ordered se- 
quences of the five occupation states, {n^,n^), reached 
as the three particles journey through the loop. The 
system transit from one occupation state to an other, 
when a particle reaches one of the two vertices of the 
loop. To make this definition clearer, we write explic- 
itly the sequences corresponding to the two scenarios, 
which chiefly rule the asymptotic dynamics. The sce- 
nario A = {(0, 1), (1,1), (1,2), (1, 1), (0,1)} is exemplified 
by the experimental pictures in Fig. |4|V. Three parti- 
cles journey simultaneously in the loop, thereby induc- 
ing a change in the particle distances. Scenario B — 
{(0,1), (1,1), (1,0), (1,1), (0,1)}, the loop is explored at 
most by two particles simultaneously. The other six traf- 
fic patterns are explicitly given, and sketched, in [14j. 
Practically, 5 is a piecewise map, which has a different 
analytical expression, Sx, for each scenario. We first 
locate the regions of the phase plane in which each sce- 
nario prevails. To do so, using Eq. [T| we compute the 
five times, t^^ , i = 1 . . . 5, at which a particle reaches 
a vertex. The linearity of the Kirchhoff laws, implies 

(i) 

that the t-^ are linear functions of Ai and A2. Con- 
sequently, the region corresponding to the scenario X 
is a polygon defined by the inequalities: t*;^-' (Ai, A2) < 
i^'''"'^'' (Ai, A2). The 8 polygons tile the phase plane as 
illustrated in Fig. [2] We can then calculate the two dis- 
tances Sx{^{n)) — (Ai(n + 1), A2(n + 1)) by computing 
the time intervals, which separate the exit of two subse- 
quent particles from the loop, and multiplying it by the 
fluid velocity outside the loops, v = + vi- Again, the 
Kirchhoff laws require the Sx to be affine functions of the 
interparticle distances: A(n-l-l) = Mx ■ X{n) + Lx, where 
the Mx and the Lx are constant matrices and constant 
vectors. Their exact but lengthy expressions are given in 
the supplemental document [14 . 

We now exploit these analytical results to give a more 
physical insight on the geometrical and dynamical prop- 
erties of the traffic dynamics. Firstly, by superimpos- 
ing the numerical trajectories on the eight regions of the 
phase plane, we notice that the asymptotic orbits are 
enclosed only in the union of the polygons A and B, 
Fig. [2j Moreover, the orbits that are enclosed in only 
one of those two regions are ellipses. To account for 
these observations, we compute the eigenvalues and the 
determinant of the Mx- Independently of the values of 
a and c, Mx is area preserving, det Mx = 1, in these 
two regions. Beyond our numerical observations, this 
central result unambiguously proves that the 3-particles 
dynamics is Hamiltonian in A and B. Furthermore, a 
tedious calculation proved that the eigenvalues of Ma 
and Mb are two complex conjugate numbers, see |14| . 
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FIG. 3: The period, r, defined from the maximum of the 
power spectra of Ai(n), is plotted versus the inverse of the 
clogging parameter c, for a — 10/9. Each point corresponds 
to a different initial condition. Full (resp. dotted) line: theo- 
retical predictions for the periods tb (resp. ta)- Inset: Phase 
portrait for c — 0.9628, the period does not depend on the 
initial condition, but the self-similar trajectories are not ellip- 
tic. 



Consequently, the orbits are necessarily self-similar el- 
lipses centered on a unique fixed point, when solely en- 
closed in A or B, in agreement with our numerical results. 
Fig. [2] In addition, the system necessarily converges to- 
ward the three Hamiltonian regions A, B, and C (re- 
gion C corresponds to the trivial case Sc = I)- Indeed, 
I detMxl takes only two different expressions elsewhere. 
I detMxl = a(l + c)/(a + c), in regions X = D,E,F,G 
and \detMx\ = (1 + c)(a + c)/[a(l + 2c)], in region 
X — H. In both cases we verify that jdetAfxl > 1, 
as 1 < a < 1 + c. This implies that, asymptotically, 
the corresponding maps yield a continuous increase of 
I All and IA2I. Therefore, as these maps are defined only 
in polygons having a finite width, we conclude that the 
system escapes from these regions as the particles flow 
through the loops. We also infer from this observation, 
that the largest invariant curve is tangent to one of the 
boundary lines of the polygon AU B, see Fig. [2] 

A second and important generic result is that the 
asymptotic traffic dynamics is time reversible. We now 
outline the demonstration of this result, which we use 
to account for the symmetry of the phase portrait with 
respect to the Ai = A2 direction. In this context, time- 
reversal corresponds to the permutation of the inter- 
particle distances: T : (Ai,A2) — (A2,Ai). Indeed, the 
last two particles that exit a loop correspond to the first 
two entering particles when reversing the flow. Saying 
that S is time-reversible thus translates into TSTS — I. 
This relation is obviously met along ID trajectories en- 
closed in only one of the two regions A or B. The cor- 
responding traffic scenarios indeed correspond to palin- 
dromic sequences of occupation states. The same re- 
sult can be also directly checked, by computing (TSx)^, 
where X = A, B, using the analytic expressions of the 
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affine maps given in [T3^. This identity is also satisfied 
for trajectories overlapping the polygons A and B as well. 
The reason for this is that T5jf (A) G X for the As be- 
longing to the invariant curves of the region X = A,B. 
The demonstration of this last result is tedious. It is 
detailed in the supplemental document [2]. In order 
to show that the global symmetry of the phase portrait 
reflects the invariance upon time-reversal symmetry, let 
us consider a ID orbit that crosses the symmetry line 
of T, at a point Ag = TAg. Noting, that = I and 
= TST, we have S-"T = TS'^T. Combining this 
relation and 5"As S"TXs yields 5"A,. = TS-"Xs. 
This last identity precisely means that the entire orbit is 
symmetric with respect to T, as any A on this orbit can 
be generated from Ag (A — S^Xg). 

We now complete this description by a brief comment 
on the properties of the mean oscillation period, t, of 
the asymptotic dynamics. For small (resp. large) c, 
the orbits are included in the region A (resp. B) only. 
Therefore, the periods ta and tb correspond to the ar- 
gument of the eigenvalues of the matrices and Mb ■ 
They are plotted versus 1/c in Fig. [3] , using the ana- 
lytic expressions given in |14| . For intermediate clog- 
ging parameters, the orbits overlap A and B. Over a 
period, ua loops are explored according to the scenario 
A and ub according to the scenario B. As the trajecto- 
ries are closed curves around the fixed point, r satisfies 
= [2Tr{nA)/TA + 27r(nB)/TB] / [(ua) + (ns)]. This 
relation implies that the oscillation period is bounded by 
TA and Tb , in agreement with Fig. [3j We also understand 
why there exists a unique period when the fixed point is 
on the boundary- line between A and B. As it includes 
the center of the ellipses, this line separates the elliptic 
orbits of both regions into two identical parts. Therefore, 
{ua) = ta/'2, and {tlb) = tb/^. This is again confirmed 
by the plots in Fig. |3]and Fig. |3] inset. 

Finally, to further confirm our theoretical predictions, 
we compare them to microfluidic experiments. Using the 
method introduced in |13j , we made a device including 20 
identical loops (L-^ — 1,675 mm, — 1.525 mm, chan- 
nel width 75 fj,ia and height 75 /im) . We monitored the 
trajectories of several triplets of identical water droplets 
advected by a continuous phase of hexadecane oil. By 
comparing the velocity of an isolated droplet in the upper 
branches and in the straight parts of the channel, we de- 
duce the experimental value of Ld — 1.2 ± 0.25 mm from 
Eq. ^ This makes possible a direct comparison between 
our experimental and our theoretical results, without any 
free fitting parameter. The evolution of Ai(n) and A2(n) 
are plotted in Fig. |4p. The grey value of each point 
codes for the traffic scenario we observed experimentally. 
Though, the fine structure of the phase portrait cannot 
be probed with a 20-loops network, an excellent agree- 
ment between our experimental and theoretical results is 
found, when considering the three generic features of the 
asymptotic-dynamics: (i) The two asymptotic-dynamics 
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FIG. 4: A: Five subsequent pictures of a typical experiment. 
Three drops (highlighted with a red dot) explore one loop 
according to the scenario A. B: Grey polygons: regions in 
which each scenario is expected from the theory. Connected 
dots: experimental trajectories recorded after 5 loops have 
been explored. The greylevel of each dot codes for the ob- 
served trafficking scenario. The blue dots correspond to a 
loop including a geometrical defect. Red line: Ai — A2. 



schemes. The distances oscillate around a fixed point 
when Ai, A2 < A„iax and the traffic scenarios are of type 
A or B only. In contrast, when Ai,A2 > Amax, we only 
observed small and non-predictible variations of the A^. 
Complete freezing was never observed due to fiuctuations 
in the droplet size, inducing differences in the droplets' 
mobility, (ii) Our model perfectly predicts the location 
of the straight boundaries between the different traffic re- 
gions, (iii) The experimental phase portrait is symmetric 
with respect to the Ai = A2 direction. 

In conclusion, combining experimental, numerical and 
theoretical tools, we have provided a comprehensive de- 
scription of the 3-body traffic dynamics. We expect that 
the generalization of our approach to coupled elementary 
traffic maps should provide a useful toolbox to design 
functional microfiuidic devices. 

We thank Laurette Tuckerman, Michael Schindler, 
Eric Lauga and Charles Baroud for valuable comments 
and discussions. We acknowledge support by C'Nano 
IdF, Sesame He de France and Paris emergence. 
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I. TRAFFIC SCENARIOS 



The eight different traffic scenarios are sketched in the figure below. 
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FIG. 1: Sketch of the traffic sequences corresponding to the eight traffic regions tlirough a single fluidic loop. 

We recall that the map 5, defined by A(n + 1) = S(\{n)) is a piecewise afiine map. Its exact expressions in the 
eight regions associated to the scenarios X = A, B , . . . , H are (setting Amax as the unit length): 

• Region A = {(Q, 1), (L 1), (1, 2), (1, 1), (0, 1)}: 



"Ai(n + 1)' 




A2(n+ 1) 





(a+c)(a+c+l) c 

a(c+l){a+2c+l) c+1 

c a(2c+l)(a+2c+l) 

c+1 (c+l)(a+c)(a + c+l) 



Ai(n) 
A2(n) 



2c-a(tt + c) + l 
c+1 

3(a + 2c+l)(a=' + (a-2)c 



(c+l){a+c)(a+c+l) 



(1) 



gionB = {(0,l),(Ll),(l,0),(l,l),(0,l)}: 



■Ai(n + iy 




A2(n + 1) 





a+c c(a + c) 

co+ti (c+l){a+c+l) 

c(a+c) a(n+2c+l)(2ca+a+2c(c+l) + l) 
(c+l)(a+c+l) (c+l)(a+c)(a+c+l)2 



AlH 



-(o-l)(o + l)2-(o-2)c2 + (-2a=' + o+2)c 
(c+l)(ti + c+l) 
a(„+2c+l)((a-l)(a+l)' + (a-2)c' + (a(2a-l)-2)c) 



(c+l)(a+c)((i+c+l) = 



(2) 



Region C = {(0, 1), (0, 0), (0, 1), (0, 0), (0, 1)}: 



"Ai(n + 1)' 




1 




■Ai(n)l 




Mn + 1) 




1 




X2{n)\ ^ 





Region D = {(0, 1), (1, 1), (1, 0), (0,0), (0, 1)}: 



Ai(n+iy 
A2(n+1) 



1 ac 

„ o(c+l) 
" o+c 



Ai(n 
A2(n) 



-a(a^+a+(a-l)c-l) + l 

o+c+1 
(a-l)a(a+c+l) 
a+c 



Region E = {(0, 1), (0, 0), (0, 1), (1, 1), (1, 0)}: 



Ai(n + 1) 
A2(n + 1) 



a+c 
c(e+l) 



(a+c)(a+c+l) 



Ai(n 

A2(n) 



(a-l)a(a+c+l) 
a+c 
:(c+l) 



(a+c)(a+c+l) 



Region F = {{0, 1), (1, 1), (1,0), (1, 1), (1, 0)}: 



■Ai(n + 1) 
A2(n + 1) 



-1 
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, o(a+2e+l) 
(a+c)(a+c+l) 



Ai(ny 
A2(n) 



a — 1 I ^ I g+l 2^ 





Region G = {(0, 1), (1, 1), (1, 2), (1, 1), (1, 0)}: 



A2(n + 1) 



_ u + c+i ac 
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Al(n) 
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Region H = {(0, 1), (1, 1), (1, 2), (0, 2), (0, 1)}: 



Ai(n + 1) 
A2(n + 1) 
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Ai(n)' 
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II. EIGENVALUES IN THE REGIONS A AND B 



The eigenvalues of the matrices Ma and Mb are complex conjugate and their modules equal to 1. The argument 
of the eigenvalues corresponds to a rotation of angle 9 in the phase plane, given by the lengthy expressions: 

• Scenario A: 



^ ^ ^ ( aia + c:)i,a + c + l)(a + 2c + 1)^4 - ^^^^^^^^^^ ^^^^^t^^^ , 
^ ^'^'^^"l c3 + (8a2 + 4a + l)c2 + 2a(a + l)(4a + l)c + 2a2(a + l)2 ' ^' 



Scenario B: 



I c3 + 2a(a+l)(3a+l)c + 2a2(a + l)2 + (2ac + c)2 I ^ ' 



We readily infer that the maximal rotation angle is bounded by 7r/2. The corresponding temporal periods are given 
by ta = 'i'T^I^A and tb = I-kIQb- 

We note also that the period of the S map takes a more simple form in the two following asymptotic limits: 
c < 1 (implying that a ~ 1), ta ~ ttv^ + (47r/y3)c-\ and c > 1, tb ~ 27r/ arctan( v'4a2 - 1) + 27r(4a^ + a - 
l)/(\/4a^ — 1 arctan(\/4a^ — l))c~^. These two asymptotic formula yield an excellent agreement with our numerical 
findings. 



III. DEMONSTRATION OF THE IDENTITY: TSTS = I 



To finalize the demonstration of this identity, we need to show that the image of the region A by the involution TS 
does not overlap the region iJ, and vice versa. To do so, let us first demonstrate that the boundary line that separates 
the regions A and B is invariant upon TS. The As that belong to this line fulfill the condition: t\ (\) = t\ {\) = 
t^'(A) = tg'(A). A straightforward calculation of these times yields the parametric equation of the boundary line: 

Al/Amax = i , , 

X /X - (a+c)(a+c+l) . I . U^i 
A2/A„,ax - - a(a+2c+l) ' + 
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Applying TSa to Eq. 11, we obtain another parametric equation for a straight hne: 




— s 



l-g 

c+1 



+ 



;(a+2c+l) 



S 



(12) 



Setting t = 



a^(a+2c+l) 



s, we readily show that the two set of equations parametrize the same line. Following the 



- (c+l){a+c+l) 



same method we show that the boundary line is invariant upon TSb as well. 

As the region A is enclosed in the triangle (q;,/3,7) defined in Fig. 2, showing that the TSa{ol = (0,0)) is located 
below the boundary line implies, that for all X € A, and lying along an invariant curve: TSTS{X) = TSaTSa{X). 
This is readily done using Eq. 1. The same reasoning apphes to region B which is enclosed in the triangle (7,5, e) 
also defined in Fig. 2. 

We then conclude that for all A e AU B and belonging to an invariant curve of S: TSTS{X) = TSaTSa(X) or 



TSTS{X) = TSbTSb{X). 



